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ABSTRACT 


An  efficient  contour-plotting  routine  is  discussed  which 
is  based  on  a  scanning  algorithm  of  Cottafava  and  LeMoli  and 
employs  bi-linear  interpolation.  An  auto-interpolation  scheme 
is  developed  which  automatically  adjusts  the  number  of  inter¬ 
polations  in  any  data  s  .are  to  produce  smooth  line  segments. 
A  program  listing  and  examples  are  given. 
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INTRODUCTION 


There  are  many  approaches  to  producing  contour  maps  on  a 
digital  computer.  Several  of  these  are  described  by  Cottafava 
and  LeMoli  (1969).  Ideally,  a  contour  routine  should  be  efficient 
on  both  computer  and  plotter,  but  most  seem  to  possess  only  one 
kind  of  efficiency.  For  quite  a  number  of  reasons,  including 
plotting  efficiency,  algorithms  of  the  "line-following"  type 
are  preferable  for  use  with  mechanical  plotters.  Cottafava  and 
LeMoli  present  a  scanning  algorithm  of  this  type  which  is  also 
very  efficient  on  the  computer. 

In  their  program  (Cottafava  and  LeMoli,  private  communica¬ 
tion)  they  assumed  a  linear  variation  between  points,  but  did 
not  do  any  interpolation  in  the  interior  of  the  data  square 
(defined  by  four  contiguous  data  points  as  vertices).  Thus  a 
plot  produced  with  their  routine  consists  entirely  of  straight 
line  segments  joined  together,  the  coarseness  depending  on  the 
spacing  of  data  points.  To  remedy  this,  the  author  developed  an 
autointerpolation  scheme  employing  bi-linear  interpolation  for 
use  in  the  interior  of  the  data  square.  With  this  scheme, 
described  in  the  present  paper,  the  number  of  interpolations 
used  in  crossing  the  data  square  is  automatically  adjusted  to 
produce  a  smooth  curve,  the  number  required  depending  on  the 
curvature  of  the  line  segment.  This  method  requires  no  additional 
storage  and  is  very  fast. 

While  the  interpolation  does  indeed  produce  smooth  line 
segments,  slope  discontinuities  sometimes  occur  on  the  edges 
(of  the  data  squares).  Tnis  is  a  limitation  of  the  bi-linear 
interpolation  law.  An  easy  and  economical  solution  to  this 
problem  is  to  obtain  a  finer  data  mesh  by  performing  a  higher- 
order  interpolation  before  entering  the  contour  routine.  This 


SCANNING  ALGORITHM 


The  scanning  algorithm  used  is  that  of  Cottafava  and 
LeMoli,  with  several  modifications  by  this  author.  A  more 
complete  discussion  than  will  be  given  here  can  be  found  in 
the  paper  cited.  For  each  contour  value,  the  procedure  is  to 
scan  the  entire  data  array  to  find  which  line  segments  are 
intersected  by  the  level  line  and  to  store  the  information  as 
flags  within  the  data  words  themselves.  A  horizontal  and 
vertical  segment  are  associated  with  each  data  point  as  shown 
in  Figure  1  (the  y-axis  is  given  its  normal  sense  here; 
Cottafava  and  LeMoli  reverse  it). 


The  flags  are  stored  as  bits  in  the  upper  part  of  the  integer 
word,  the  data  being  normalized  to  positive  integers  by  a  linear 
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transformation.  (In  deference  to  readers  with  machines  which 
cannot  perform  masks,  the  use  of  masks  in  the  programs  has 
been  limited  to  an  inessential  role,  where  they  can  easily  be 
replaced  by  a  test  and  subtract.  Where  masking  statements  are 
available,  the  flags  can  be  conveniently  stored  as  the  least- 
significant  bits  of  the  floating-point  mantissa  of  each  data 
word.  This  results  in  a  considerable  simplification  of  the 
p  rog  ram. ) 

The  contours  are  traced  in  a  second  scanning  operation, 
each  flag  being  erased  as  it  is  found.  This  procedure  makes  it 
easy  to  find  all  the  branches  of  the  contour  level  and  is  in 
large  part  responsible  for  the  efficiency  of  the  program.  Each 
line  is  traced  square  by  square  by  a  local  scan  which  checks 
all  the  edges  of  the  data  square  in  fixed  order  (counter¬ 
clockwise  beginning  with  the  right  edge)  to  find  the  contin¬ 
uation  of  the  line.  This  procedure  runs  into  trouble  only  in 
the  case  of  an  interior  saddle  point  (square  crossed  twice  by 
the  same  contour),  and  in  this  case  there  is  an  easy  solution 
based  on  the  interpolation  method.  j 
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INTERPOLATION 


In  order  to  define  the  behavior  of  the  contour  line 
inside  each  data  square,  we  must  make  art  assumption  about  the 
behavior  of  the  function  inside  the  square.  Lacking  any  special 
information  in  the  general  case,  we  assume  bi-linear  variation 
as  the  simplest  general  variation.  (Using  a  higher  order  inter¬ 
polation  here  would  also  cause  serious  difficulties  with  the 
scanning  procedure.)  That  is,  we  assume 


F(x,y)  =  A  +  0x  +  ay  +  6xy 


(1) 


referred  to  a  local  coordinate  system  with  origin  at  A,  as 
shown  in  Figure  2. 

y 


t  v 

(0,1)  C  S - "v - ^  o 
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Equation  (1)  is  the  Lagrange  2x2  interpolation  formula  and 
is  also  equivalent  to  a  Taylor's  series  expansion,  with  the 
assumption  of  constant  first  derivatives  on  each  edge.  Making 
our  square  of  unit  dimensions,  the  derivatives  are 


C-A 


i  :  B-A 
Y  -  D-C 


(2) 


a  a  Y-S 


With  these  definitions,  it  is  easy  to  see  that  (1)  reduces  to 
the  correct  values  at  the  corners  and  reduces  to  ordinary 
linear  interpolation  on  each  edge.  Our  contour  segment  is 
therefore  the  locus  F(x,y)  =  v,  the  value  of  the  contour. 

From  (1)  we  obtain  either 


v-A-Sx 
a+ ox 


(3) 


o  r 


y  =  v-h-ay 
3  +6y 


(4) 


with  0«x<l  and  0_<y£l  .  These  expressions  are  easily  computed. 

In  practice,  we  choose  the  number  of  interpolations,  subdivide 
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Ax  or 

Ay, 

and 

use 

convenient 

to 

use 

edge , 

and 

(4) 

for 

either  (3)  or  (4),  respectively.  It  is 
(3)  for  a  line  terminating  on  a  vertical 
a  line  terminating  on  a  horizontal  edge. 
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AUTO- INTERPOLATION 


We  approximate  our  ideally  smooth  curve  by  a  series  of 
chords.  If  we  choose  the  number  of  chords  in  each  square  so 
that  the  maximum  deviation  from  the  ideal  curve  is  on  the 
order  of  the  basic  plotter  increment,  then  we  obtain  as  smooth 
a  curve  as  we  can  with  no  wasted  time.  We  take  as  our  auto- 
interpolation  criterion,  the  maximum  pe rpend i cu 1  a r  distance 
from  the  curve  to  the  straight  line  between  the  end  points. 

The  perpendicular  distance  from  a  point  to  a  line  is 


• (x)  =  (y 


m  x  -  b  )  /"I 


(5) 


from  elementary  geometry.  We  consider  the  case  of  a  segment 
terminating  on  the  left  edge  of  the  square,  and  for  this  case 
there  are  just  two  distinct  possibilities,  as  shown  in  Figures 
3  and  4.  All  other  possibilities  can  be  found  by  reflections 
and  a  rotation. 
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For  case  I  the  slope  is 


m 


1 


v-C  .  v-C 
C-A  '  C-D 


-Y/'* 


For  case  1 1  , 


v-A-jj  _  v-A  _  _  i  (  v  -  A  )  +  i 
m2  u+o  c*  a ( a*  h ) 


For  both  cases  the  intercept  is  b  =  (v-A)/i.  The  "brute  force" 
calculation  of  the  maximum  of  (5)  is  messy,  but  a  transformation 
simplifies  the  algebra.  Define 


n 


1 


i(v-C) 

ay 


n2  -  6 /a, 

n  =  —  Ax  where  Ax 

a 


((  v  —  C ) / Y  ,  case  I 
J  ,  case  II 


(6) 


Then 


i(v-A) 


e  =  Ytl  +  ',] 


a 


+  3 


6  [  1  +  ^-]  + 


which  leads  to  the  relation 


m2  ( 1 ^  ; 2  ) 


( 1 +n1 ) 


(7) 


It  can  also  be  shown  that 


y-b 


m  1  ( 1  +  • ,  ]  ) 
1  +  o„x 


x 


and 


y 


iZ.  = 

dx 


( 1 +n]  ) 
( l+n2x)2 


The  location  of  the  maximum  of  (5)  is  given  by  the  condition 


i 

i 


c  '  (  x  )  -  0  =  y  1  (  x  )  -  m 


(hats  will  be  used  to  refer  to  the  maximum),  i.e.. 


m  i  (  1  +  'i  i  ) 
(!+•  2*)2 


(8) 
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Thus 


y-b  =  m ( 1 +  r\^x ) * 


and 


Vm^+1  e  =  y  -  b  -  mx 


mn2x 


or 


e 


(9) 


where  x  is  found  from  (8).  For  case  I,  m  =  m^  ,  and 


■  '-2  h/T^T  -  n 

For  case  II,  m  =  m  2 ,  and  using  (7)  in  (8)  gives 

x2  =  n^  tVl+r| 2  "  1  ^ 

which  has  the  same  form.  Putting  these  results  in  (9)  gives 
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m 


for  both  cases.  This  is  conveniently  written  in  the  form 


tpr  (VT^7  -1)  )  do) 

Since  (10)  is  symmetric  in  x  and  y,  it  can  be  easily  seen  to 
apply  to  the  case  of  a  segment  terminating  on  the  bottom  edge 
of  the  square  also,  provided  that  we  define  n  :  <5Ay/S  in  this 
case.  Thus  all  possible  cases  are  contained  in  (10).  A  bound 
can  be  placed  on  c  by  considering  n-*00  and  Ax-*-Ay-*-l  ,  namely 
~sl/*T,  which  agrees  with  geometrical  intuition. 

Knowing  how  to  calculate  the  maximum  deviation  of  the  curve 
from  the  straight  line  between  endpoints,  we  use  this  to  estimate 
the  number  of  segments  needed  to  approximate  the  curve  to  any 
desired  accuracy.  To  do  this,  we  consider  the  case  of  a  circular 
arc,  Figure  5. 

L  j 

I 

I  r 

I 

1 

Figure  5.  Deviation  from  a  chord 
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f 


2 

It  is  easily  shown  that  d=L  /8r,  provided  L/r<<l,  the  important 
point  being  the  proportionality  d  L^.  Since  L  “  1/N,  approximately, 
where  N  is  the  number  of  segments  needed  in  that  square,  we  take 


N  -  1  +  /|e|/d 


where  d  equals 
increment,  d  = 
scheme  appears 


the  allowable  deviation.  For  plotters  with  a  c 
.001"  is  satisfactory.  This  auto- i nterpol ation 
to  work  quite  well. 


5  mil 
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SINGULARITIES  AND  S ADDL E - PO I  NTS 


Evidently  something  peculiar  happens  with  (10)  if  o+l<0. 
This  can  be  seen  to  be  connected  with  a  singularity  in  y  or  x, 
equation  (3)  or  (4).  The  existence  of  a  singularity  in  y(x) 
within  the  interval  0£x<_l  is  implied  by  the  condition  that  (D-B) 
and  a  have  opposite  signs,  and  the  existence  of  a  singularity  in 
x(y)  within  the  interval  0_<y_<1  is  implied  by  the  condition  that 
y  and  £  have  opposite  signs.  If  both  singularities  exist,  then 
the  square  has  an  interior  saddle  point  located  at  the  inter¬ 
section.  An  example  is  shown  in  Figure  6.  A  very  convenient 
criterion  for  making  connections  in  a  saddle  square  is  that 
contours  shoul d  never  cross  a  singularity.  The  condition 
n+l<0  can  be  easily  shown  to  imply  a  wrong  connnection  in  a 
saddle  square.  By  far  the  easiest  solution  to  this  problem  is 
just  to  check  for  a  wrong  connection  and  to  resume  scanning 
the  square  if  it  exists.  At  most  three  tries  will  be  necessary 
to  make  the  correct  connection,  and  since  saddle-points  should 
be  relatively  rare,  this  is  a  small  price  to  pay  for  such  a 
simple  procedure  that  guarantees  correct  connections. 
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RECTANGULAR  MESH 


Should  it  be  desired  to  plot  data  cells  as  rectangles 
instead  of  squares,  this  is  easily  done  by  stretching  one  axis 
in  the  calls  to  PLOT.  Defining  r  as  the  ratio  of  x  to  y  scale 
factors,  i.e.  the  rectangle  has  length  r  in  the  x-direction 
and  1  in  the  y-direction,  elementary  trigonometry  gives  for  the 
modi f i ed  deviation 


V 


l  +  r2m2 
1+m2 
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EXAMPLES 


Two  simple  examples  of  plots  produced  by  the  routine 
are  shown  in  figures  7  and  8.  Figure  7  was  produced  from  real, 
deterministic,  data,  and  Figure  8  from  random  numbers.  Data 
points  are  marked  by  ticks  along  the  borders  of  the  plots.  Each 
plot  is  based  on  only  24  data  points,  and  the  large  number  of 
slope  discontinuities  indicates  the  need  for  a  more  refined 
data  mesh. 
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APPENDIX 
PROGRAM  LISTING 

Written  in  FORTRAN-63,  a  programming  language  of  the 
CDC  1604  computer.  Integer  words  assumed  at  least  33  bits  long. 


2-0 


UJ  07  /I 


c 

L 

c 

t 

l 

c 

c 

c 

c 

c 

c 


c 

c 

c 


1 

V 

3 

!> 


7 

8 


18 


30 

35 

<0 

50 
6o 
1 0 
80 
100 


lUPHOU  I’„fc  CONIOS*  R.N.F.  8.C.  NU.NUF. IV. 1A,  Mfc  UlMT.UEV  ) 

•  I  RE  R  5  I  0  "  r  « ft  > »  F|N1)»H|,  iV(NU2»M)»  I  A  1  R  ) 

•  *TA  UhDK  ■  lU9aOOOOOOOfl’'  <,V6H  "  tQOUUUOOUOU®  1 

•LOT  CONTOURS  F!i»J>  ■  FU.X)  •  0*  J*1»N  »  J*i.M 

IV  »  1A  ARE  AUXILIARY  ARRAYS.  ray  eUUlVAtbNLb  !F.IV>» 

■UST  RAVE  ND2  t 

•  OF V*  is  AUTOMNTERPOLATION  PARAHfclbK  *  AHHHUX.  alLO«ABlE  LfcV. 

•  iNCHbS  )  .  SUREST  DEV  ■  .uoi 

•  LGORIThm  BY  CUTTAFAvA  f  lFPOOT.  POL  I TELM  CO  U1  R1LACO. 

1RANSCR I RED  \  REVISED  BY  N*RK  RIRTH.  JANUARY  xV7i> 
interpolation  algbrithr  SY  RA«N  *1H1r 

•ORNALIZF  DATA 

•HAX  ■  F  S  FRI*  •  V. 

■  0  5  «  ■  1 #  H 
10  i  I  i  1,N 

IF <  FU.jl.GT.FMAX  1  l»2 
•MAX  «  FU.Jl  •  GO  10  5 

IF  C  FU.jl.LT.FNlN  >  3.5 
•NIN  «  F!l.J> 

QONT I NUg 

•  ■  lt*e/<FMAX*FMi N»  »  8  ■  •  ►  M I N»A  ».s 

•R  ■  Pi  •  R  *  l  »  Ni  «  N  ♦  1  »  NRb  «  N  •  2 

■0  8  w  ■  1 ■ H 
10  7  1  ■  l.N 

I  V ( 1 ♦ 1 , JR )  ■  A»FU,JE-1)  *  8 
IV(I.WR)  ■  IV<NP2*JR>  ■  o 

<R  »  jR  -  . 

•  ALL  ERASE!  nP*.  I  V  > 

•  ALL  ERASE!  nP7» I  V  1 1, R.J )  i 

•  AC  ■  HEI(.HT/<*«i*  •  DEVS  *  DEV  /  FAC 

•ALL  FACTOR!  FAC  i 

LOOP  CVga  CONTUUR  LfcVELS 
10  50U  L  *  1.8 

•  0  ■  l VOL  ■  A *L ( L A  ♦  B 
•relirin»rt  scan,  set 

|A!j9  •  Mftl  -  I  VOL 

■  '»•>• . 

OONTINUE 

•0  loc  I  *  2.Nl 

II  »  I  ♦  1 

•o  loo  J  *  Z.M1 
IF!  I.EQ.N1  >  Eq.Aq 
1C  •  IVllI.j)  »  I  VOL 
IF!  IC  )  40,35 
1C  ■  i  « 

IF!  XSlGNf !*, i A ! JI)»IC 
1 V  !  I  «  w  1  ■  IVll»J>  ♦  IVER 
IF!  J.bQ.Mj  )  l0n.70 

IF!  XblGNPll. JA< >  eo»8u»loo 
IVII.v)  «  Iv!l«J>  ♦  IHOR 
IAIJI  ■  IC 


FLa<*s 

IV«^»J>  •  IV‘*»J>  *  x 

IVUJ.J  I  *  IVUl.Jl  ♦  1 
)  5u«5g.6(j 


-o 


/ 


0  HICK  l.p  UF  CONTOU"  Lt»tl 

QO  2nu  .)  *  ?»" 

L  BCTTOI*  ETufc 

I  1 v  «  1 Y <2 . j ) . A\n. .NrT . j vg* 

lf(  Ilv.'.e.inUH  >  i2u» 1 3ft 
l2o  OALL  o"A‘Lm<  Ng? .  1  v  >  *0,uEYi>  ) 

l.  TOP  FU*E 

lJl)  If  <  I  *  t  Nj  •  J ) .  ufc ,  I  hPH  )  i4g»20Q 
1 4  u  QAlL  '>»l»J*l*2»  Nn2» IV. wU.OtVb  ) 

?00  0 Q *J T  I  N  U  g 
C  HlGMT  tQP.fc 

*0  igu  I  *  2»1'1 

1F(  fVlt.Hj^i.St.lVgH  )  23„»j00 

?5U  OALL  oHANLM<  1, Hi. 2*3.  N02  »  1  V  ,  VO,  Ufe  *  S  ) 
3  j  ij  OONTINUg 

L  LEFT  tUOE  »  INlgSiOK 

D°  5oq  J  ■  2.H 

•»  PI  »  c  /  J 

10  50U  I  «  2,N 

JF(  )  5an.J50 

J5B  BALL  b  H  4  nL  M  (  ItJ.i.KH.  ND2  *  I  V  .  V  0  »  U  t  V  S>  ) 

5  o  0  COMlNUg 

OALL  PACTURI  X.  I 
■ETurn 
•  NO 


5tJ  r  -  I  I  !t  iKA.*„H<  1-L^(  LISP.A*  Nu*  i  V  ,  Vlj,  Jt  V  S  ) 

i  1  hcnS  I  0‘  1V<N0*1> 

JATA  llHOWa  1  .  LOOUri  0  0  0UB  1  *  1  *  ^tKI  2u0  UO  0  0  VUUl)b  1  *  1  1  VM*  3o  0  U  0  U  u  U  0  0  u  ° 
0  POLIO"  O'.fc  RRANCh  OF  C11NT0UK  L  t  V  t  L 

0 

|PFN  *  3  *  IP  •  !  1  |HtM  A  JP  *  J  *  jREM 

IP  *  V  -  2  A  TP  a  J  -  2  A  it  *  U 

NJSP  *  LtSP  *  ASSIGN  i  ru  JAjL  A  GO  Tu  200 

1  IPFN  *  2 

SO  TO  <5.4,3).  K 

3jS*vsJ«i  A  GO  10  S 

4  IS  X  |  3  i  .  1 

2  IPt  I  V < I . j»l ) .oe. 1 Vt«  1  10*12 

10  4  •  J*1  A  MlSP  «  2  A  ASSIGN  &  lo  jail  S  GO  !0  2UU 

1 2  IbV  «  I V< !♦,  ,J».AND.  .NOT. I VfcH 
IP  C  IoV.Gfc.lHOK  )  2U.25 

20  1  *  1*1  A  K I AP  ■  1  »  ASSIGN  S  U  J*iL  A  go  10  2  Q  U 

22  IF<  I  * < I ,J> .Gfc. IVfcS  I  30*35 

3U  *ISP  *  2  »  ASSIGN  J  to  JAIL  A  GO  1u  20U 

32  IF  <  I V< I , J) .Gt. Ih8h  )  4fl. 9U 

4  o  4 1 SP  «  i  A  ASSIGN  <t  TO  JAJL  A  GO  10  20u 

*0  I F  <  K  )  100*1*0 

100  RETURN 
C 

o  HEMOvt  flags,  interpolate,  anj  ploi  poini 

190  I  »  INfcM  »  J  ■  JPE"  a  IVtl.JJ  ■  JV(1»J)»IVEH 

PISP  ■  2  A  ASSIGN  mo  10  JAIL 

2oO  *0  *  v  •  2  A  TO  ■  I  -  2 

40  *  JP  A  10  a  IP 

IF (  J.fcQ.JP  )  410*240 
210  40  *  vP  ♦  1 
220  IF<  I*fc0.lP  )  43n*24u 
23o  10  »  IP  *  l 

2  4  Q  4  a  IA  ■  1 v< I • J) ,ANU. .NOT, I VH 
a  *  lb  a  l  V(  I  *OQ)  *ANI1.  .NOT*  I  VK 
0  «  IL  a  1VUU»J)*ANO..NOT;IVH 
I  a  It  a  IV(  [U*  JCl!  .AND.  .NO'  .  IVK 

4LF  S  C  -  A  »  UfcT  a  B  -  A  A  OfcU  ■  U-C  *  BfcT 

80  TO  *220.22^*  "ISP 

25o  *F  *  <VO-A*/pEI  *  XO  A  0*  «  *F  «  XP  $  IF  »  TO 

SIGN  ■  1>J  -  I  A  UT  a  X  a  SlUN*<YP-Y0> 

I TA  a  Dy»UEL  /  BET  A  GO  10  2/u 

26e  IF  ■  ‘VO-AI/ALF  ♦  YO  a  DY  a  IF  •  YP  S  XF  «  XO 

SIGN  »  JO  -  J  A  Q  A  a  A  a  SiGN»«Xp-X0> 

•TA  a  OX*CEL  /  ALF 
270  IF (  ETA  )  2an.40r 

Jflfl  ITAP  «  ETA  *  i* 

C  SADOLfc-POlNT  TfcST 

I F (  ETAp  )  290*292.300 

?9o  I  »  IS  *  J  •  JS  »  II  *  H  ♦  j 

SO  TO  <i5*25.32),  IT 
292  IPS  •  i .  A  GO  To  Sin 

300  IPS  ■  DX«DYa(SURTMETAP)-l« »**4  ‘  < b | A*SOHTf I UX aux*D Y »D V > > 

BIO  II  ■  1.5  ♦  SQH1F(ABSF!EPS>/UEVS) 

INI  a  N l  * , 999V  »  IT  i  , 

BO  TO  <320*35,.  >.  "ISP 


U  3  02  /I 


J2u  OY  «  T  /  FN[ 

S  »  V  L  •  A  -  A  t  ►  *  V  »  UP  «  *AL ►  *U  T 

9  *  Htl  •  DEL*t  »  t)L>  »  PfcL*UT 

J3q  »  =  Y  -  nr 

If  (  Y  )  4  0  0,4  'U,  34(1 

540  F  *  F  -  "F  »  G  *  G  •  UG  A  X  ■  MG 

CALL  PLOT!  xn»X,YU*S!GPj.Y,  IPfcN  ) 

30  TO  3  ? o 
35u  I*  s  *  /  Fnj 

P  s  VL-A  *  hfT*»  *  UP  ■  *dfcl*UX 

a  =  ALA  *  Dr  L  *  *  »  DG  ■  W£L»UX 

J40  I  »  X  -  n* 

IF  (  X  )  4  0  0 , 4  ,U.  J/Q 

B7Qp3F-0F  »  G  ■  0,  -  UG  I  Y  ■  F/G 

8ALL  FL0T(  XO»l>lGl»***  T0*Y,  IP£N  ) 

90  TO  JAU 

4  Q  Q  Iw(l.w)  •  -  [H0fi«MlaP 

IS  *  IP  ■  I  »  JS  •  JP  «  J  i  XP  *  XP  i  YP  *  Tf 

«ALU  PLOT  t  XP»TP,  IPEN  ) 

90  TO  JAIL,  <1»3M»5»ioo* 

INC 


